### File    : corstange-bjmes-replication.r
##  Author  : Daniel Corstange (daniel.corstange@columbia.edu)
##  Created : Mon 27 Nov 2017 --- 11:43 AM
##  Notes   : Replication file for "The Syrian Conflict and Public
##            Opinion Among Syrians in Lebanon"
##            (British Journal of Middle Eastern Studies)

### preliminaries: load libraries and data

library(dplyr)
library(lattice)
library(ggplot2)
library(psy)
library(psych)

Syr <- read.csv("corstange-replication-syria.tsv", sep = "\t")
AB3 <- read.csv("corstange-replication-arab-barometer.tsv", sep = "\t")

### plots
#### factions

## gather the data and fix the factor ordering
Dat_faction <-
    with(Syr,
         data_frame(support = as.vector(prop.table(table(faction_all))),
                    group   = as.factor(levels(faction_all))))
Dat_faction$group <-
    with(Dat_faction,
         factor(group,
                levels(group)[rev(c(3, 4, 1, 2, 5))]))

## plot it
plot_faction <- 
    dotplot(group ~ support,
            data   = Dat_faction,
            dat    = Dat_faction,
            xlim   = c(-.02, .62),
            xlab   = "Percentage",
            scales = list(x = list(cex    = 0.75,
                                   at     = seq(0, 1, .10),
                                   labels = seq(0, 100, 10))),
            panel  = function(x, y, dat, ...) {
                ## leaders
                panel.segments(x0  = rep(-1000, 5),
                               x1  = as.numeric(x),
                               y0  = as.numeric(y),
                               y1  = as.numeric(y),
                               lty = 1,
                               col = "grey50")
                ## Islamist extra segments
                panel.segments(x0  = dat$support[grep("Islamist",
                                                      dat$group)],
                               x1  = sum(dat$support[grep("Islamist",
                                                          dat$group)]),
                               y0  = c(3, 2),
                               ## y0  = c(2, 3),
                               y1  = 2.5,
                               lty = 3,
                               col = "grey67")
                ## opposition extra segments
                panel.segments(x0  = c(dat$support[grep("National",
                                                        dat$group)],
                                       sum(dat$support[grep("Islamist",
                                                            dat$group)])),
                               x1  = sum(dat$support[grep("Islamist|National",
                                                          dat$group)]),
                               y0  = c(4, 2.5),
                               y1  = 3.25,
                               lty = 3,
                               col = "grey67")
                ## extra dots     
                panel.xyplot(x   = c(sum(dat$support[grep("Islamist",
                                                          dat$group)]),
                                     sum(dat$support[grep("Islamist|National",
                                                          dat$group)])),
                             y   = c(2.5, 3.25),
                             pch = 19,
                             cex = 1.50,
                             col = "grey67")
                ## extra values
                panel.text(x      = c(sum(dat$support[grep("Islamist",
                                                           dat$group)]),
                                      sum(dat$support[grep("Islamist|National",
                                                           dat$group)])),
                           y    = c(2.5, 3.25) + .25,
                           labels = round(100 * c(sum(dat$support[grep("Islamist",
                                                                       dat$group)]),
                                                  sum(dat$support[grep("Islamist|National",
                                                                       dat$group)])), 0), 
                           cex    = .50,
                           col    = "grey67")
                ## extra dot labels
                panel.text(x      = c(sum(dat$support[grep("Islamist",
                                                           dat$group)]),
                                      sum(dat$support[grep("Islamist|National",
                                                           dat$group)])),
                           y    = c(2.5, 3.25) - .25,
                           labels = c("Islamists", "Opposition"),
                           cex    = .75,
                           adj    = 0,
                           col    = "grey67")
                ## dots
                panel.xyplot(x   = as.numeric(x),
                             y   = as.numeric(y),
                             pch = 19,
                             cex = 1.50,
                             col = "black")
                ## values
                panel.text(x      = as.numeric(x),
                           y      = as.numeric(y) + .25,
                           labels = round(as.numeric(x) * 100, 0),
                           cex    = .50,
                           col    = "grey50")
            })

## plot to pdf
pdf(file    = "plot-faction.pdf",
    width   = 4 * 1.618,
    height  = 4,
    family  = "Palatino",
    onefile = FALSE)
plot_faction
dev.off()

#### room density

## gather data
Dat_room <-
    with(Syr,
         data_frame(dens    = c(syr_room_density,
                                leb_room_density),
                    group   = rep(factor(ifelse(faction_none == 1, NA,
                                         ifelse(faction_govt == 1, 2, 1)),
                                         labels = c("Opposition",
                                                    "Government")),
                                  2),
                    country = rep(factor(1:2,
                                         labels = c("Syria",
                                                    "Lebanon")),
                                  each = length(syr_room_density))))

## plot it
plot_room <-
    ggplot(subset(Dat_room, !is.na(group)),
           aes(x     = group,
               y     = dens,
               color = group,
               fill  = group)) + 
    ## draw this first to get the caps on the whiskers
    stat_boxplot(geom = "errorbar",
                 width = 0.25) +
    ## draw this second to get a "filled" box and no outliers
    stat_boxplot(fill          = "white",
                 notch         = TRUE,
                 notchwidth    = 0.75,
                 outlier.shape = NA,
                 varwidth      = TRUE) +
    ## underplot the actual points
    geom_jitter(width  = 0.25,
                height = 0.10,
                size   = 0.75,
                alpha  = 0.075) +
    ## overlay the actual box and whiskers
    geom_boxplot(alpha         = .25,
                 notch         = TRUE,
                 notchwidth    = 0.75,
                 outlier.shape = NA,
                 varwidth      = TRUE) +
    ## tinker with scales
    scale_color_manual(values = c("blue3", "red3")) +
    scale_fill_manual(values = c("blue3", "red3")) +
    labs(x = NULL,
         y = "People per Bedroom") +
    coord_flip(ylim = c(0,8)) +
    facet_wrap(~ country, nrow = 2) +
    ## tinker with themes
    theme_bw() +
    theme(legend.position = "none")

## plot to pdf
pdf(file    = "plot-room-density.pdf",
    height  = 4,
    width   = 4 * 1.618,
    family  = "Palatino",
    onefile = FALSE)
plot_room
dev.off()

#### age

## gather data
Age <-
    data.frame(vals  = c(as.vector(prop.table(table(Syr$age_bins_05))),
                         as.vector(aggregate(faction_all ~ age_bins_05,
                                             FUN  = function(x) prop.table(table(x)),
                                             data = Syr)$faction_all)) * 100,
               bins  = rep(as.factor(levels(Syr$age_bins_05)), 6),
               group = rep(factor(1:6,
                                  labels = c("Overall",
                                             "None",
                                             "Government",
                                             "Nationalists",
                                             "Domestic\nIslamists",
                                             "Foreign\nIslamists")),
                           each = nlevels(Syr$age_bins_05)))

## plot it
plot_age <-
    ggplot(subset(Age, group == "Overall" & bins != "60-65"),
           aes(x = bins, y = vals)) +
    ## background barchart for overall age distribution
    geom_bar(stat  = "identity",
             color = "grey85",
             fill  = "grey95") +
    ## foreground lines for factions
    geom_line(data    = subset(Age, (!(group %in% c("Overall", "None")) &
                                     bins != "60-65")),
              mapping = aes(x     = as.numeric(bins),
                            y     = vals,
                            size  = group,
                            alpha = group,
                            color = group)) +
    ## annotate the lines
    annotate(geom  = "text",
             x     = c(3, 1.5, 8, 6),
             y     = c(44, 32, 11, 29),
             label = levels(Age$group)[!(levels(Age$group) %in%
                                         c("Overall", "None"))],
             size  = c(2.5, 3, 2.5, 3),
             alpha = c(.50, 1, .50, 1),
             color = c("grey33", "black", "grey33", "black")) +
    ## tinker with scales
    labs(x = "Age Cohort",
         y = "Percent of Age Cohort") +
    scale_y_continuous(limits = c(0, 50)) +
    scale_size_manual(values = c(1.75, 1, 1.75, 1)) +
    scale_alpha_manual(values = c(1, .33, 1, .33)) +
    scale_color_manual(values = c("black",
                                  "grey33",
                                  "black",
                                  "grey33")) +
    ## tinker with themes
    theme_bw() +
    theme(legend.position = "none")

## plot to pdf
pdf(file    = "plot-age.pdf",
    height  = 5,
    width   = 5,
    family  = "Palatino",
    onefile = FALSE)
plot_age
dev.off()

#### sex

## gather data
Sex <-
    data.frame(vals  = as.vector(aggregate(faction_all ~ female,
                                           FUN = function(x) prop.table(table(x)),
                                           data = subset(Syr,
                                                         faction_all != "No Preference"))[[2]][,-5]) * 100,
               sex   = factor(1:2,
                              labels = c("M", "F")),
               group = rep(factor(1:4,
                                  labels = c("Domestic Islamists",
                                             "Foreign Islamists",
                                             "Government",
                                             "Nationalists")),
                           each = 2))
Sex$focus <-
    with(Sex,
         ifelse(group %in% c("Nationalists",
                             "Foreign Islamists"),
                TRUE, FALSE))
Sex$group <-
    with(Sex,
         reorder(group, -vals))
Sex$xpos <-
    with(Sex,
         as.numeric(group) +
         .125 * ifelse(sex == "M", -1, 1))

## plot it
plot_sex <-
    ggplot(Sex,
           aes(x      = xpos,
               y      = vals,
               label  = sex,
               color  = focus)) +
    geom_segment(aes(x    = xpos,
                     xend = xpos,
                     y    = 0,
                     yend = vals)) +
    geom_point(size  = ifelse(Sex$focus, 8, 6),
               shape = 21,
               fill  = "white") +
    geom_text(size = ifelse(Sex$focus, 4, 3)) +
    labs(x = NULL, y = "Percent of Sex") +
    scale_x_continuous(labels = levels(Sex$group),
                       minor_breaks = NULL) +
    scale_color_manual(values = c("grey50", "black")) +
    coord_cartesian(ylim = c(7.5, 42.5)) +
    theme_bw() +
    theme(legend.position = "none",
          axis.text.y = element_text(size  = 7),
          axis.text.x = element_text(size  = 7,
                                     color = c("grey50",
                                               "black",
                                               "black",
                                               "grey50")))

## plot to pdf
pdf(file    = "plot-sex.pdf",
    height  = 5,
    width   = 5,
    family  = "Palatino",
    onefile = FALSE)
plot_sex
dev.off()

#### religion

## gather data
Dat_rel <-
    data_frame(group  = c("Domestic Islamists",
                          "Foreign Islamists",
                          "Government",
                          "Nationalists"),
               piety  = aggregate(rel_index ~ faction_all,
                                  FUN  = function(x) prop.table(table(x)),
                                  data = subset(Syr,
                                                faction_all != "No Preference"))[[2]][,"2"],
               relpol = aggregate(relpub_index ~ faction_all,
                                  FUN  = function(x) prop.table(table(x)),
                                  data = subset(Syr,
                                                faction_all != "No Preference"))[[2]][,"3"])
Dat_rel$group <-
    with(Dat_rel, reorder(group, -relpol))

## plot it
plot_rel <- 
    dotplot(group ~ relpol,
            xlim   = c(.175, .625),
            xlab   = "Percentage",
            data   = Dat_rel,
            dat0   = Dat_rel,
            key    = list(space        = "top",
                          columns      = 2,
                          padding.text = 0,
                          points       = list(pch  = c(21, 21),
                                              col  = c("black",
                                                       "black"),
                                              fill = c("black",
                                                       "grey95"),
                                              cex  = 1.00),
                          text         = list(labels = c("Public Life",
                                                         "Private Life"),
                                              cex    = 0.75)),
            scales = list(x = list(cex    = 0.75,
                                   at     = seq(0, 1, .10),
                                   labels = seq(0, 100, 10)),
                          y = list(cex = 0.75,
                                   col  = c("grey50",
                                            "grey50",
                                            "grey50",
                                            "black"))),
            panel  = function(x, y, dat0, ...) {
                panel.segments(x0  = rep(-1000, 8),
                               x1  = c(dat0$relpol, dat0$piety),
                               y0  = rep(as.numeric(y), 2) +
                                   .15 * c(rep(1, 4), rep(-1, 4)),
                               y1  = rep(as.numeric(y), 2) +
                                   .15 * c(rep(1, 4), rep(-1, 4)),
                               lty = 1,
                               col = rep(c("grey75",
                                           "grey75",
                                           "grey50",
                                           "grey75"), 2))
                                           ## "grey50"), 2))
                panel.xyplot(x    = c(dat0$relpol, dat0$piety),
                             y    = rep(as.numeric(y), 2) +
                                 .15 * c(rep(1, 4), rep(-1, 4)),
                             pch  = 21,
                             cex  = rep(c(1, 1, 1.25, 1), 4),
                             col  = rep(c("grey50",
                                          "grey50",
                                          "black",
                                          "grey50"),
                                        2),
                             fill = c("grey50",
                                      "grey50",
                                      "black",
                                      "grey50",
                                      "grey95",
                                      "grey95",
                                      "white",
                                      "grey95"))
            })

## plot to pdf
pdf(file    = "plot-rel.pdf",
    width   = 4 * 1.618,
    height  = 4,
    family  = "Palatino",
    onefile = FALSE)
plot_rel
dev.off()

#### attentiveness

## gather data
Dat_engaged <- 
    data_frame(group  = c("Domestic Islamists",
                          "Foreign Islamists",
                          "Government",
                          "Nationalists"),
               eng    = aggregate(engaged_index ~ faction_all,
                                  FUN  = function(x) prop.table(table(x)),
                                  data = subset(Syr,
                                                faction_all != "No Preference"))[[2]][,"3"])
Dat_engaged$group <-
    with(Dat_engaged, reorder(group, eng))

## plot it
plot_engaged <-
    dotplot(group ~ eng,
            data   = Dat_engaged,
            dat    = Dat_engaged,
            xlim   = c(-.025, .525),
            xlab   = "Percentage",
            scales = list(x = list(cex    = 0.75,
                                   at     = seq(0, 1, .10),
                                   labels = seq(0, 100, 10)),
                          y = list(cex    = 0.75,
                                   col    = c("black",
                                              "grey50",
                                              "grey50",
                                              "black"))),
            panel  = function(x, y, dat, ...) {
                ## leaders
                panel.segments(x0  = rep(-1000, 4),
                               x1  = as.numeric(x),
                               y0  = rep(as.numeric(y), 4),
                               y1  = rep(as.numeric(y), 4),
                               lty = 1,
                               col = c("grey75",
                                       "grey50",
                                       "grey75",
                                       "grey50"))
                ## dots
                panel.xyplot(x   = as.numeric(x),
                             y   = as.numeric(y),
                             pch = 19,
                             cex = c(1, 1.50, 1, 1.50),
                             col = c("grey50",
                                     "black",
                                     "grey50",
                                     "black"))
                ## values
                panel.text(x      = as.numeric(x),
                           y      = as.numeric(y) + .20,
                           labels = round(as.numeric(x) * 100, 0),
                           cex    = .50,
                           col    = "grey50")
            })

## plot to pdf
pdf(file    = "plot-engaged.pdf",
    width   = 4 * 1.618,
    height  = 4,
    family  = "Palatino",
    onefile = FALSE)
plot_engaged
dev.off()

#### age appendix

## gather data from UN-sponsored Vulnerability Assessment
UN_age <-
    data.frame(bins   = factor(1:15,
                               labels = c("0-4",
                                          "5-9",
                                          "10-14",
                                          "15-19",
                                          "20-24",
                                          "25-29",
                                          "30-34",
                                          "35-39",
                                          "40-44",
                                          "45-49",
                                          "50-54",
                                          "55-59",
                                          "60-64",
                                          "65-69",
                                          "70+")),
               male   = c(20, 18, 13, 8, 5,
                          6,  9,  7,  4, 3,
                          2,  2,  1,  1, 1) / 100,
               female = c(17, 16, 12, 9, 9,
                          10, 8,  6,  4, 3,
                          2,  2,  1,  1, 1) / 100)
UN_age <-
    with(UN_age,
         data.frame(bins = bins,
                    age  = c(male, female),
                    sex  = rep(factor(1:2,
                                      labels = c("Male",
                                                 "Female")),
                               each = nrow(UN_age))))

## gather survey data
Syr$agebins5 <-                         # bin survey data
    with(Syr,
         cut(age,
             c(-Inf, seq(24, 64, 5), Inf),
             labels = c("20-24",
                        "25-29",
                        "30-34",
                        "35-39",
                        "40-44",
                        "45-49",
                        "50-54",
                        "55-59",
                        "60-64",
                        "65-69")))
Syr$agebins5 <-                         # bin survey data to match UN data
    Syr$agebins5_alt <-
        with(Syr,
             cut(age,
                 c(-Inf, seq(4, 69, 5), Inf),
                 labels = c("0-4",
                            "5-9",
                            "10-14",
                            "15-19",
                            "20-24",
                            "25-29",
                            "30-34",
                            "35-39",
                            "40-44",
                            "45-49",
                            "50-54",
                            "55-59",
                            "60-64",
                            "65-69",
                            "70+")))

## get the survey data into dataframes
Age <-
    data.frame(bins = levels(Syr$agebins5),
               age  = c(with(subset(Syr, female == 0),
                             prop.table(table(agebins5))),
                        with(subset(Syr, female == 1),
                             prop.table(table(agebins5)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Syr$agebins5))))
Age_alt <-
    data.frame(bins = ordered(1:length(levels(Syr$agebins5_alt)),
                              labels = levels(Syr$agebins5_alt)),
               age  = c(with(subset(Syr, female == 0),
                             prop.table(table(agebins5_alt))),
                        with(subset(Syr, female == 1),
                             prop.table(table(agebins5_alt)))),
               sex  = rep(factor(1:2,
                                 labels = c("Male", "Female")),
                          each = length(levels(Syr$agebins5_alt))))

## put UN and survey data together (adults only)
UN_adult <-
    subset(UN_age,
           !(bins %in% c("0-4",
                         "5-9",
                         "10-14",
                         "15-19")))
UN_adult$age <- 
    c(with(subset(UN_adult, sex == "Male"),
           age / sum(age)),
      with(subset(UN_adult, sex == "Female"),
           age / sum(age)))
Survey_adult <- 
    subset(Age_alt,
           !(bins %in% c("0-4",
                         "5-9",
                         "10-14",
                         "15-19")))
Adult_combo <- rbind(UN_adult, Survey_adult)
Adult_combo$src <-
    rep(factor(1:2,
               labels = c("UN Vulnerability Assessment", "Survey")),
        each = nrow(Adult_combo) / 2)
Adult_combo$age <-
    with(Adult_combo,
         ifelse(sex == "Male", age, -age))

## plot it
plot_age_both <-
    ggplot(data = Adult_combo,
           aes(x    = bins,
               y    = age,
               fill = sex)) +
    geom_hline(yintercept = 0) +
    geom_bar(stat = "identity",
             data = subset(Adult_combo, sex == "Male"),
             alpha = 0.50) +
    geom_bar(stat = "identity",
             data = subset(Adult_combo, sex == "Female"),
             alpha = 0.50) +
    annotate(geom  = "rect",
             xmin  = (nrow(Adult_combo) / 4) - 0.33,
             xmax  = (nrow(Adult_combo) / 4) + 0.33,
             ymin  = 0.125 * c(-1, 1) - 0.05,
             ymax  = 0.125 * c(-1, 1) + 0.05,
             color = "grey50",
             fill  = "white") +
    annotate(geom     = "text",
             x        = nrow(Adult_combo) / 4,
             y        = -0.125,
             fontface = "italic",
             color    = "grey50",
             label    = "Female") +
    annotate(geom     = "text",
             x        = nrow(Adult_combo) / 4,
             y        = 0.125,
             fontface = "italic",
             color    = "grey50",
             label    = "Male") +
    scale_fill_manual(values = c("darkred", "darkblue")) +
    scale_y_continuous(limits = c(-0.25, 0.25),
                       breaks = seq(-0.3, 0.3, 0.1),
                       labels = abs(round(seq(-0.3, 0.3, 0.1), 1))) +
    labs(x    = "Age",
         y    = "Proportion") +
    coord_flip() +
    facet_wrap(~ src) +
    theme_bw() +
    theme(legend.position = "none")

## plot to pdf 
pdf(file    = "plot-age-appendix.pdf",
    height = 4,
    width  = 8,
    family = "Palatino",
    onefile = FALSE)
plot_age_both
dev.off()

#### province appendix

## gather data: survey and UNHCR
## UNHCR figures compiled from:
##   UNHCR Lebanon - Beirut Country Office. Syria Refugee Response
##   Lebanon: Places of Origin of Syrian Refugees Registered in
##   Lebanon. 31 March 2015.
Province <-
    data.frame(province = levels(factor(Syr$province)),
               survey_n = as.vector(table(factor(Syr$province))),
               unhcr_n  = c("Al-Hasaka"    = 2711 + 15233 + 12351 + 618,
                            "Al-Raqqa"     = 5678 + 43411 + 15850,
                            "Al-Suwayda"   = 251 + 666 + 124,
                            "Daraa"        = 23386 + 30613 + 28354,
                            "Deir al-Zour" = 17169 + 3395 + 5094,
                            "Dimashq"      = 55105,
                            "Halab"        = 4618 + 23882 + 26423 + 28898 + 11140 + 10608 + 133659 + 3094,
                            "Hama"         = 16048 + 50585 + 7681 + 11229 + 956,
                            "Homs"         = 1089 + 3072 + 151189 + 15150 + 55705 + 19940,
                            "Idlib"        = 4412 + 67666 +10462 + 24061 + 45247,
                            "Latakia"      = 3464 + 919 + 40 + 520,
                            "Quneitra"     = 10512 + 31,
                            "Rif Dimashq"  = 19175 + 62755 + 14502 + 25305 + 12524 + 2679 + 11440 + 2326,
                            "Tartus"       = 1749 + 25 + 49 + 1181 + 345))
Province$survey_prop <-
    with(Province,
         survey_n / sum(survey_n))
Province$unhcr_prop <-
    with(Province,
         unhcr_n / sum(unhcr_n))
Province$survey_vs_unhcr <-
    with(Province,
         survey_prop - unhcr_prop)
Province$province <-
    with(Province,
         reorder(province, -survey_prop))

## plot it
plot_province <-
    ggplot(data = Province,
           aes(x = province,
               y = survey_prop)) +
    geom_bar(stat = "identity",
             aes(x = province,
                 y = survey_prop),
             fill  = "darkblue",
             alpha = .50) +
    geom_bar(stat = "identity",
             aes(x = province,
                 y = -unhcr_prop),
             fill  = "darkred",
             alpha = .50) +
    geom_hline(yintercept = 0) +
    geom_line(aes(x = as.numeric(province),
                  y = survey_vs_unhcr),
              color = "darkorchid4") +
    geom_point(aes(x = as.numeric(province),
                   y = survey_vs_unhcr),
               color = "darkorchid4") +
    geom_label(x        = length(Province$province),
               y        = -.125,
               size     = 3,
               fontface = "italic",
               color    = "grey50",
               label    = "UNHCR Estimate") +
    geom_label(x        = length(Province$province),
               y        = .125,
               size     = 3,
               fontface = "italic",
               color    = "grey50",
               label    = "Survey Sample") +
    scale_y_continuous(limits = c(-0.25, 0.25),
                       breaks = seq(-0.3, 0.3, 0.1),
                       labels = abs(round(seq(-0.3, 0.3, 0.1), 1))) +
    labs(x    = "Province",
         y    = "Proportion") +
    coord_flip() +
    theme_bw()

## plot to pdf
pdf(file    = "plot-province-appendix.pdf",
    height = 4,
    width  = 4,
    family = "Palatino",
    onefile = FALSE)
plot_province
dev.off()

### tables and in-text stats
#### Arsal bounds

## combine domestic and foreign Islamists
Syr$faction_arsal <-
    with(Syr,
         factor(ifelse(faction_all == "No Preference", 1,
                ifelse(faction_all %in% c("Foreign Islamists",
                                          "Domestic Islamists"), 2,
                ifelse(faction_all == "Nationalists", 3,
                ifelse(faction_all == "Government", 4, NA)))),
                labels = c("No Preference",
                           "Islamists",
                           "Nationalists",
                           "Government")))

## actual proportions in the full sample
prop.table(table(Syr$faction_arsal))

## extreme bounds: what if all Baalbeki govt supporters were Islamists?
baalbekis <-
    with(subset(Syr, district == "Baalbek"),
         table(faction_arsal))
baalbekis_extreme <- baalbekis
baalbekis_extreme[c("Government", "Islamists")] <-
    c(0, sum(baalbekis[c("Government", "Islamists")]))
## ... and here are the proportions under these extreme conditions
prop.table(with(subset(Syr, district != "Baalbek"),
                table(faction_arsal)) +
           baalbekis_extreme)

## less extreme: replace Baalbekis with proportions of rest of sample
prop.table(with(subset(Syr, district != "Baalbek"),
                table(faction_arsal)) +
           round(with(subset(Syr, district != "Baalbek"),
                      prop.table(table(faction_arsal))) *
                 sum(Syr$district == "Baalbek"), 0))

#### respondent type

## what type of respondent:
## - fully randomized
## - a male who replaced a female because SHE refused to participate
## - a male who replaced a female because HE refused to allow her to participate

prop.table(table(Syr$respondent_type))

#### room density

## room density in Syria
with(Syr,
     quantile(syr_room_density,
              na.rm = TRUE))
with(subset(Syr, faction_all == "Government"),
     quantile(syr_room_density,
              na.rm = TRUE))
with(subset(Syr, !(faction_all %in% c("Government",
                                      "No Preference"))),
     quantile(syr_room_density,
              na.rm = TRUE))

## room density in Lebanon
with(Syr,
     quantile(leb_room_density,
              na.rm = TRUE))
with(subset(Syr, faction_all == "Government"),
     quantile(leb_room_density,
              na.rm = TRUE))
with(subset(Syr, !(faction_all %in% c("Government",
                                      "No Preference"))),
     quantile(leb_room_density,
              na.rm = TRUE))

## deterioration of room density from Syria to Lebanon
with(Syr,
     prop.table(table(log(leb_room_density /
                          syr_room_density) > 0)))
with(subset(Syr, faction_all == "Government"),
     prop.table(table(log(leb_room_density /
                          syr_room_density) > 0)))
with(subset(Syr, !(faction_all %in% c("Government",
                                      "No Preference"))),
     prop.table(table(log(leb_room_density /
                          syr_room_density) > 0)))

## difference in room density variance between Syria and Lebanon
with(Syr,
     var.test(log(leb_room_density),
              log(syr_room_density)))
with(subset(Syr, faction_all == "Government"),
     var.test(log(leb_room_density),
              log(syr_room_density)))
with(subset(Syr, !(faction_all %in% c("Government",
                                      "No Preference"))),
     var.test(log(leb_room_density),
              log(syr_room_density)))

## summary statistic table details
## - columns are 1Q, median, 3Q, geometric mean, geometric variance
## - rows are full sample, govt, oppn, nationalists, Islamists,
##   Domestic Islamists, Foreign Islamists
round(rbind(with(Syr,                   # Syria room density
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, faction_all == "Government"),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, !(faction_all %in% c("Government",
                                                  "No Preference"))),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, faction_all == "Nationalists"),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, faction_all %in% c("Foreign Islamists",
                                                "Domestic Islamists")),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, faction_all == "Domestic Islamists"),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density)))),
            with(subset(Syr, faction_all == "Foreign Islamists"),
                 c(quantile(syr_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(syr_room_density), na.rm = TRUE)),
                   exp(var(log(syr_room_density), na.rm = TRUE)),
                   sum(!is.na(syr_room_density))))),
      digits = 2)

round(rbind(with(Syr,                   # Lebanon room density
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, faction_all == "Government"),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, !(faction_all %in% c("Government",
                                                  "No Preference"))),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, faction_all == "Nationalists"),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, faction_all %in% c("Foreign Islamists",
                                                "Domestic Islamists")),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, faction_all == "Domestic Islamists"),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density)))),
            with(subset(Syr, faction_all == "Foreign Islamists"),
                 c(quantile(leb_room_density,
                            probs = c(.25, .50, .75),
                            na.rm = TRUE),
                   exp(mean(log(leb_room_density), na.rm = TRUE)),
                   exp(var(log(leb_room_density), na.rm = TRUE)),
                   sum(!is.na(leb_room_density))))),
      digits = 2)

#### education

## illiteracy rates by sex
with(subset(Syr, female == 0),          # men
     prop.table(table(educ_illit)))
with(subset(Syr, female == 1),          # women
     prop.table(table(educ_illit)))

## table for education outcomes
educ_table <- 
    round(rbind(with(Syr,
                     c(prop.table(table(educ)))),
                with(subset(Syr, faction_all == "Government"),
                     c(prop.table(table(educ)))),
                with(subset(Syr, !(faction_all %in% c("Government",
                                                      "No Preference"))),
                     c(prop.table(table(educ)))),
                with(subset(Syr, faction_all == "Nationalists"),
                     c(prop.table(table(educ)))),
                with(subset(Syr, faction_all %in% c("Domestic Islamists",
                                                    "Foreign Islamists")),
                     c(prop.table(table(educ)))),
                with(subset(Syr, faction_all == "Domestic Islamists"),
                     c(prop.table(table(educ)))),
                with(subset(Syr, faction_all == "Foreign Islamists"),
                     c(prop.table(table(educ)))),
                with(AB3,
                     prop.table(table(educ)))),
          digits = 2)
colnames(educ_table) <-
    c("L", "M", "H")
rownames(educ_table) <-
    c("Full Sample",
      "Government",
      "Opposition",
      "Nationalists",
      "Islamists",
      "Domestic Islamists",
      "Foreign Islamists",
      "Arab Barometer")
educ_table

#### sect

## Sunni Arab vs everyone else
with(Syr,
     prop.table(table(minority)))

## Kurdish speakers
with(Syr,
     prop.table(table(kurdish_speak)))

## sects
with(Syr,
     prop.table(table(sect)))

## table
with(subset(Syr, faction_all != "No Preference"),
     table(minority,
           faction_all == "Government"))
with(subset(Syr, faction_all != "No Preference"),
     prop.table(table(minority,
                      faction_all == "Government")))

#### religion

## table for religious outcomes
rel_table <- 
    round(rbind(with(Syr,
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, faction_all == "Government"),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, !(faction_all %in% c("Government",
                                                      "No Preference"))),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, faction_all == "Nationalists"),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, faction_all %in% c("Domestic Islamists",
                                                    "Foreign Islamists")),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, faction_all == "Domestic Islamists"),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(subset(Syr, faction_all == "Foreign Islamists"),
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)), 
                       prop.table(table(rel_index)),
                       prop.table(table(relpub_family_bi)),
                       prop.table(table(relpub_index)))),
                with(AB3,
                     c(prop.table(table(pray_daily)),
                       prop.table(table(quran_weekly)),
                       prop.table(table(rel_index)),
                       rep(NA, times = 5)))),
          digits = 2)
colnames(rel_table) <-
    c("L", "H", "L", "H", "L", "M", "H", "L", "H", "L", "M", "H")
rownames(rel_table) <-
    c("Full Sample",
      "Government",
      "Opposition",
      "Nationalists",
      "Islamists",
      "Domestic Islamists",
      "Foreign Islamists",
      "Arab Barometer")
rel_table

## chi-squared test for opposition groups: high vs not-high piety
with(Syr,
     chisq.test(table(rel_index == 2,   # high vs not
                      faction_all)[,c("Nationalists",
                                      "Domestic Islamists",
                                      "Foreign Islamists")]))

## correlation of public and private religion
with(Syr,
     cor(cbind(rel_index,
               relpub_index),
         use    = "pairwise.complete.obs",
         method = "spearman"))


## Cronbach's alpha
with(Syr,                               # all 4: .83
     cronbach(cbind(relpub_econ,
                    relpub_pol,
                    relpub_crime,
                    relpub_family)))

with(Syr,                               # no family: .92
     cronbach(cbind(relpub_econ,
                    relpub_pol,
                    relpub_crime)))

## correlations
relpub_cor4 <-                             # all 4: mean .65
    with(Syr,
         cor(cbind(relpub_econ,
                   relpub_pol,
                   relpub_crime,
                   relpub_family),
             use    = "pairwise.complete.obs",
             method = "spearman"))
relpub_cor3 <-                             # no family: mean .88
    with(Syr,
         cor(cbind(relpub_econ,
                   relpub_pol,
                   relpub_crime),
             use    = "pairwise.complete.obs",
             method = "spearman"))
relpub_cor4_tri <- relpub_cor4[lower.tri(relpub_cor4)]
fisherz2r(mean(fisherz(relpub_cor4_tri)))
range(relpub_cor4_tri)
relpub_cor3_tri <- relpub_cor3[lower.tri(relpub_cor3)]
fisherz2r(mean(fisherz(relpub_cor3_tri)))
range(relpub_cor3_tri)

#### attentiveness

## table for attentiveness outcomes
attentive_table <- 
    round(rbind(with(Syr,
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, faction_all == "Government"),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, !(faction_all %in% c("Government",
                                                      "No Preference"))),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, faction_all == "Nationalists"),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, faction_all %in% c("Domestic Islamists",
                                                    "Foreign Islamists")),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, faction_all == "Domestic Islamists"),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(subset(Syr, faction_all == "Foreign Islamists"),
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index)))),
                with(AB3,
                     c(prop.table(table(understand)),
                       prop.table(table(interest)),
                       prop.table(table(know)),
                       prop.table(table(engaged_index))))),
          digits = 2)
colnames(attentive_table) <-
    c("L", "H", "L", "M", "H", "L", "M", "H", "L", "M", "H")
rownames(attentive_table) <-
    c("Full Sample",
      "Government",
      "Opposition",
      "Nationalists",
      "Islamists",
      "Domestic Islamists",
      "Foreign Islamists",
      "Arab Barometer")
attentive_table

## chi-squared tests: govt vs opposition
with(subset(Syr, faction_all != "No Preference"),
     chisq.test(table(engaged_index,
                      faction_all == "Government")))
with(subset(Syr, faction_all != "No Preference"),
     chisq.test(table(understand,
                      faction_all == "Government")))
with(subset(Syr, faction_all != "No Preference"),
     chisq.test(table(interest,
                      faction_all == "Government")))
with(subset(Syr, faction_all != "No Preference"),
     chisq.test(table(know,
                      faction_all == "Government")))

## Cronbach's alpha
with(Syr,                               # .82
     cronbach(cbind(understand,
                    interest,
                    know)))

## correlations (Spearman's rho)
attentive_cor <-                        # mean = .66
    with(Syr,
         cor(cbind(understand,
                   interest,
                   know),
             use    = "pairwise.complete.obs",
             method = "spearman"))
attentive_cor_tri <-
    attentive_cor[lower.tri(attentive_cor)]
fisherz2r(mean(fisherz(attentive_cor_tri)))
range(attentive_cor_tri)

## correlations (polychoric)
attentive_polychoric <-                 # mean = .85
    with(Syr,
         polychoric(cbind(understand,
                          interest,
                          know)))$rho
attentive_polychoric_tri <-
    attentive_polychoric[lower.tri(attentive_polychoric)]
fisherz2r(mean(fisherz(attentive_polychoric_tri)))
range(attentive_polychoric_tri)

### (Emacs local variables)

## Local Variables:
## eval: (outline-minor-mode 1)
## eval: (font-lock-mode 1)
## outline-regexp: "^\###+ "
## outline-level: (lambda ()
##                  (cond
##                   ((looking-at "^### ")     1)
##                   ((looking-at "^#### ")    2)
##                   ((looking-at "^##### ")   3)
##                   ((looking-at "^###### ")  4)
##                   ((looking-at "^####### ") 5)
##                   (t 1000)))
## eval: (font-lock-add-keywords
##        nil
##        '(("^### .*" 0 (aref outline-font-lock-faces 0) t)
##          ("^#### .*" 0 (aref outline-font-lock-faces 1) t)
##          ("^##### .*" 0 (aref outline-font-lock-faces 2) t)
##          ("^###### .*" 0 (aref outline-font-lock-faces 3) t)
##          ("^####### .*" 0 (aref outline-font-lock-faces 4) t)))
## eval: (outline-hide-sublevels 1)
## End:
